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The recently introduced H-function extension of the Negative Binomial Distribution is investi- 
gated. The analytic form of P„ is rederived by means of the Mellin transform. Applications of the 
HNBD are provided using experimental data for P„ in e^e~, e^p, inelastic pp and non-diffractive 
pp reactions. 



, In multiparticle phenomenology we celebrate the 25th anniversary year of the invention of two ideas which turned 
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out to be highly influential: the Koba-Nielsen-Olesen (KNO) scaling of the multiplicity distributions P„ [1] and their 
Negative Binomial Distribution (NBD) type parametrization [2]. Both subjects received a great deal of attention 
during the past 25 years, see refs. [3-5] for review articles. 

A popular way of introducing the NBD as a model of P„ is the Poisson transform defined by 
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Pn= / i^{z) — e-'-dz. (1) 
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I For the NBD the continuous probability density i}}{z) is the gamma distribution providing the asymptotic KNO scaling 
form of Pn [3]. At high energies the agreement between the NBD and observations is often not quite satisfactory 
therefore it is of interest to look for possible generalizations. We shall consider here an extension of the NBD that 
04 ' can be obtained via the scaling transformation z — s- z^/^^ of the gamma variate with scaling exponent /i > 0. The 
modified ^(z) can be written in the form 

l> , 

O . ^i^,) ^ A^'^z^'=-iexp(-[Az]'^) (2) 



which is the generalized gamma distribution [6] with shape parameter A; > and scale parameter A > 0. Its moments 
are given by 



For = 1 the ordinary gamma distribution is recovered and the Poisson transform of Eq. (2) yields the NBD, 



where {n) is the average multiplicity. For /i 7^ 1 the basic properties of the Poisson transformed generalized gamma 
density have been determined in two previous Letters [7,8]. Since the analytic form of P„ can be expressed in terms 
of the H-function of Fox [9] (see the Appendix for a summary) we will call the distribution as HNBD for short. The 
main goal of the present paper is the analysis of various multiplicity data which are known to be in disagreement with 
the ^ = 1 special case. 

In ref. [7] the analytic form of the HNBD was obtained utilizing the Laplace transform of the Fox function given 
by the equation pair (A. 6-7). This is a straightforward method since n!P„ in Eq. (1) is the Laplace transform of the 
product z^^{z) which can be expressed in terms of H(a;) for the generalized gamma distribution. To illustrate the 
simplicity of manipulating H-functions, first we provide a slightly different approach to the analytic form of the newly 
developed probability law. The Laplace and Mellin transforms of a function f{x), x G (0, 00), are defined by 







C{f{x),r) = / e-^-f{x)dx (5) 

and by 

M{f{x),s) = / x'-'f{x)dx (6) 
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respectively. They are related to each other through 



M{C {f{x),r), s} = T{s) M{f{x), 1 - s} 
which can be justified by Fubini's theorem [10]. Prom the moments of tp{z) given by Eq. (3) one obtains 

r{k+{s + n-l)/n) 1 



M{z"^{z),s} = 
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Xs+n-l 



for the Mellin transform of z""tjj{z) and the use of Eq. (7) yields 
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M{jC (^Xz),r),s} 



A"r(fc) 



r{s)r{k + n/ii- s/iJ.) X". 



(7) 
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Comparison with the MeUin transform (A. 8) of the Fox function shows that C {z"'tp{z),r = 1), providing n!P„, is 
expressible in terms of H(a;) with the following parameters: 
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A"r(A;) 
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(1 - fc - n//z, l//z) 
(0,1) 



for > 1 



(10) 



where the restriction on /i follows from the existence condition «; > of H(a;) discussed in the Appendix. Thus, using 
identity (A. 4), one arrives at 



Pn 
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(1 - k, 1//Z) 

(n,l) 



for > 1 
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and with the help of identity (A. 3) we can write the analytic form of P„ in terms of a legal H-function also for /n < 1: 



Pr,.= 
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for < ;U < 1. 



(12) 



A reparametrization of Eqs. (11) and (12) using (n) in place of the scale parameter A is more convenient. The factorial 
moments of P„ are equivalent to the ordinary moments of il){z) for Eq. (1) therefore (n) = (z) and Eq. (3) yields 



A = 



r(fc + i/M) 
r(fc)(n) • 



(13) 



For unit scaling exponent fi the HNBD reduces to the pure NBD given by Eq. (4) which is the fi = 1 marginal case of 
Eq. (11) for (n) < k and of Eq. (12) for (n) > k. This can be deduced by comparing Eq. (10) with (A. 9) and recalling 
that iFo(a, — ; x) yields the binomial function (1 — x)~°' for \x\ < 1. 

At first glance it may seem strange that a three-parameter discrete probability law requires three different expressions 
for Pn each having a restricted domain of validity over the parameter space. But this complication is a reasonable 
price for the high degree of flexibility of the HNBD. It involves as a special or limiting case the Poisson transform of 
many widely used probability densities such as the 

• gamma distribution for = 1 and A; > 0, 

• WeibuU distribution for fc = 1 and /U > 0, 

• chi distribution for = 2 and fc > 0, 

• Pareto distribution for fc ^ and — > oo, 

• log-normal distribution for — > and fc — > oo 

to mention but a few. Therefore one may think the HNBD as a family of different discrete probability laws rather 
than a single distribution. Since the infinite divisibility of ijj{z) is preserved by P„ for Eq. (1) we can use Bondesson's 
theorem [11] to deduce that the HNBD is infinitely divisible if < /z < 1, otherwise this feature of P„ is violated. 

Before turning our attention to data analysis let us consider briefly the computational aspects of the HNBD. 
Unfortunately, the Fox function is not available in software packages for the computation of special functions. A 
straightforward way to calculate P„ in fitting procedures is the evaluation of the integral Eq. (1) numerically with 
f{x) given by Eq. (2), requireing (z) = 1 and changing z to {n)z in the Poisson weight. For — > and fc — > oo the 
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HNBD converges to the Poisson transformed log- normal distribution which lacks a representation in terms of known 
functions. To include this important limiting case into the fits Eq. (2) should be written in the following form [7]; 

^i^) = -pji^^^ exp[p-^to-e"'] forp^O (14) 

where w = p{\nz + a)/a with new parameters a > 0, p and a. This reparametrization of the generalized gamma 
density allows /U < as well (/x = p/a) and the log-normal distribution can be mapped to the origin in p, 



tp{z) = — exp 
v27r az 



{In z-af 



2a2 



for p = 0. (15) 



Restricting the location parameter a to 



rinr(p- + a/p)-lnr(p-) for p ^ 

1^ -cr^/2 for p = ^ ' 

Pn is computed by inserting Eqs. (14) and (15) into Eq. (1), changing z to {n)z in the Poisson weight and evaluating 
the integral by the Romberg method or Gaussian quadrature [12] or using the numerical integration package of 
Mathematica [13]. 

Applying the HNBD to experimental data we shall consider typical reactions and energies where the pure NED is 
claimed to be in contradiction with observations. Significant deviations from = 1 are expected in 

1. inelastic pp reactions up to ISR energies, 

2. deep-inelastic e+p scattering at HERA energies, 

3. e+e~ annihilations at LEP energies, 

4. non-diffractive pp collisions at top SPS energy ^/s = 900 GeV. 

Items 2 and 3 have already been investigated in refs. [7,8] and here only a few additions will be made. But let us first 
consider inelastic pp reactions, perhaps the earliest typical example for the failure of NBD fits. Without the aim of 
completeness we have analysed multiplicity distributions in the cm. energy range ^/s = 10.7 — 546 GeV [14] including 
the inelastic pp data of the UA5 Collaboration. In the fitting procedure numerical evaluation of the integral Eq. (1) 
was carried out as described above. After some initial playing with the parameters it was foimd that the k = \ special 
case of the HNBD, i.e. the Poisson transformed WeibuU distribution, provides good description of the investigated 
data. The results of fits are collected in Table 1. It is worth noticing that the best-fit value of ji slowly decreases 
with increasing energy, further, /U > 1 for each data set. According to some preliminary results the k = 1 special case 
produces similar quality fits for inelastic Tr^p and K'^p reactions. 

The Poisson transformed WeibuU law works successfully in deep- inelastic e^p scattering at HERA energies as well. 
In ref. [8] the midtiplicity distributions measured by the HI Collaboration [15] were compared to the WeibuU case of 
the HNBD neglcicting systematic errors of P„ . To avoid possible misunderstanding we repeated the analysis with the 
inclusion of both sources of experimental uncertainties. The average multiplicity (n) was fixed at its observed value 
and only ji was treated as free parameter in the fitting procedure. The results are collected in Table 2, the quoted 
errors of fj, are statistical and systematic. The success of the k = 1 special case of the HNBD (even with omitted 
systematic errors of P„) suggests that the inelastic pp and deep-inelastic e^p multiplicity distributions are similarly 
shaped. This seems to be confirmed by the observation [15] that the HI data and the best-fit NBD deviate from each 
other at small multiplicities. Due to the diffractive fraction of events, the same happens for inelastic pp data too. In 
both reactions the rise of P„ at small n is less rapid than that of the NBD and this property is well described by the 
Poisson transform of the WeibuU distribution with /x > 1. 

A completely different type of deviation from /i = 1 arises in non-diffractive pp collisions at = 900 GeV. The 
multiplicity distributions measured by the UA5 Collaboration [16] are at present the highest energy published data 
for Pn. The Negative Binomial fits are satisfactory in narrow pseudorapidity intervals but for \ri\ > 2.5 the NBD 
fails to reproduce the observed shape of P„. According to widespread opinion the failure of the NBD is caused by 
the much quoted shoulder structure. In narrow bins this is not visible but for \ri\ > 2.5 some kind of structure can 
indeed be seen in the tail of P„. However, the dominant source of discrepancy between the NBD and observations 
has nothing to do with the shoulder. It was demonstrated already in [16] that the real difficulty for wide 77-intervals 
is the pronounced narrow peak of the heavy tailed distributions. Performing HNBD fits reveals clearly that the /x = 1 
special case is unable to reproduce a highly skewed shape which is so characteristic of the UA5 data for \r]\ > 2.5. But 
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letting /i to vary freely in the fits it turns out that the ^ 0, fc — > oo limiting case of the HNBD, i.e. the Poisson 
transformed log-normal distribution, yields reasonable description of P„ possessing a narrow peak and extended tail. 
The rcsiilts for full phase-space and wide 77- windows are collected in Table 3. Although systematic deviations between 
the experimental and theoretical P„ inevitably occur due to the shoulder structure (the HNBD is unimodal) the overall 
shape of the multiplicity distributions can be reproduced remarkably well by the two-parameter Poisson transformed 
log-normal law. For the full phase-space data the fit is illustrated in Fig. 1. 

Interestingly, the log-normal limit of the HNBD was found to be successful also in e+e~ annihilations at the 
peak. In ref. [7] the results of three-parameter HNBD fits were presented favouring /ii « for the full phase-space 
multiplicity distributions. Here we provide the outcome of two-parameter Poisson transformed log-normal fits, sec 
Table 4. The analysed data sets now include the Aleph data with updated systematic errors and the Opal data 
corresponding to s/s = 161 GeV [17]. As is seen the quality of fits are satisfactory for each experiment. Similar 
results can be obtained for the Aleph data measured in central rapidity windows [18], these are quoted in Table 5. We 
have to mention that the Delphi data in y-intervals exhibit significant shoulder structure which can not be reproduced 
by unimodal distributions and therefore the HNBD fails as well. 

The importance of the log-normal distribution to approximate P„ in e+e" and pp collisions has already been 
stressed in the literature [19]. One may think that the observed success of the Poisson transformed log-normal law 
is only another manifestation of the above cited results. But the widely known and accepted evidences for the log- 
normality of Pn at present energies are not rarely questionable. Their comparison to our findings (which indicate that 
log-normality emerges asymptotically) will be treated in a separate paper [20] . 

In conclusion, we have investigated the H-function extension of the NBD obtained by the Poisson transform of the 
generalized gamma distribution Eq. (2). The analytic form of P„ was rederived via Mellin transform. Since Eq. (2) 
involves as special and limiting cases many classical probability densities, the HNBD is expected to obey high degree 
of flexibility. Our applications confirmed this expectation. Fitting the HNBD to various multiplicity data claimed to 
be in disagreement with the pure NBD we have obtained satisfactory results. According to these findings, two major 
types of departure from fi=l arise in multiparticlc production. For inelastic pp and deep-inelastic e^p scattering the 
dominant source of discrepancy is the less steep rise of Pn at small multiplicities caused by the diffractive component 
of the underlying dynamics (at least for pp data). The multiplicity distributions can be reproduced successfully by 
the WeibuU case of the HNBD with > 1. The other p ^ 1 type behaviour arises in non-diffractive pp collisions at 
\fs = 900 GeV where the NBD is incompatible with the highly skewed shape of Pn for I77I > 2.5. These heavy tailed 
distributions with a narrow peak are best described if /z ^ and k ^ 00, i.e. by the log-normal limit of the HNBD. 
Despite of the absence of skew multiplicity curves, the same type of departure from the pure NBD has been observed 
also in e+e" annihilations at LEP energies. This remarkable ubiquity of the /x — > limit for the e+e~ and pp data at 
top energies has a natural explanation based on renormalization group arguments for asymptotic KNO scaling [20]. 
It will be interesting to see how precisely can the Poisson transformed log-normal law reproduce the forthcoming pp 
data at y/s = 1800 GeV at Tevatron. 

The observed success of the HNBD indicates that the majority of existing multiplicity data can be interpreted on 
the basis of the asymptotic KNO scaling form of the Negative Binomial Distribution with the extension of its validity 
to positive powers of the scaling variable. It is hoped that in the next 25 years the HNBD will provide a helpful 
unifying framework to study the NB regularity and the origin of possible deviations. 
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Appendix: The H-function of Fox 

The Fox function Hp '"(a;) is defined by the Mellin-Barnes type integral 



iim.n 

Hp,q 



Hp,q 



(ai,ai), • • • , (ap,ap) 
(6i,/3i),...,(6g,/3,) 



A{s)B{s) 
2m Jc C{s)D{s) 



x" ds 



with X ^0 and 



^(s) = n"=ir(l-a,+a,s), 



= Wj^m+i r(i - &i + D{s) = u'j^n+i n^j - 



(A.1) 



(A.2) 



Above m,n,p,q are integers satisfying < n < p and I < m < q. In the cases n = 0, q — m and p — n (A.2) has 
to be interpreted as i?(s) = 1, C(s) = 1 and D{s) = 1, respectively. The parameters (ai, . . . , Op) and (6i, ■ . ■ ,bq) are 
complex, whereas parameters (ai, . . . , ap) and . . . , are positive numbers. They are restricted by the condition 



that aj{hh + i') ^ PiMj - 1 - A) for z/, A = 0, 1, 



1, . . . , m; J = 1, . . . , n. The contour C in the complex s plane 



is such that the points s = {bh + i')/ Ph and s = {aj — 1 — u)/ aj lie to the right and left of C respectively while C 
extends from s = oo — ik to s = oo + ik where fc is a constant with k > |Im bh\/Ph- 

The Fox function makes sense only if the following two existence conditions are satisfied: 

1. X ^ and k > with k = Ylj=i ~ X]j=i 

2. K = and < \x\ < 1/p with p = n^^i a"' nj=i PJ^' ■ 

Under these conditions H(a;) is an analytic function for a; ^ 0, in general multivalued, one- valued on the Riemann 
surface of In x. 

Elementary properties of the Fox function very useful to manipulate them are 



iim.n 



p, q 



(ap,ap) 

{P'Pl ^p) 

{b„l3,) 



un.m 
•^q.P 



'^P.q 



(1 - b„ Pg) 

(1 - Up, Up) 



{ap + cap, Up) 

{bq + Cf3g, Pq) 



(A.3) 
(A.4) 



Ljni.n 
^P,q 



-Hp.q 



(ttp, ap/c) 
ib„ /3Jc) 



c> 0. 



(A.5) 



Many of the integral transforms of H (x) yield Fox functions again with altered parameters. For example, the Laplace 
transform is given by the equation pair 



c{H-'"^{cx),r} = -H;+^:ri 



(1 



(0,1), (l-Op-ap, ap) 



and 



m,n+l 



£{H-;^(ca;),r}=-H^T 



(1, 1), (ap + Up, ap) 

{bq+Pq, Pq) 



for < K < 1 



for K > 1 



with scale parameter c > 0. The Mellin transform of H(x) reads 

nr=i nbj + M uu r(i - - 



M{H'^:"^{cx),s} 



UUm+i r(i - bj - M uU+i + 



(A.6) 



(A.7) 



(A.8) 



An important special case of H(a;) is the generalized hypergeometric function pFq(a;). Its relation to the Fox function 
is given by 



n^^ir(«, 



tti , . . . , ap 
bi,...,bq 



"p,q+l 



UPil 

"q+l,P 



(1,1), ib„l) 
{ap, 1) 



(1 - ap, 1) 
(0,1), {l-bq,l) 



foY p < q, ov p = q + 1 and \x\ < 1. 



(A.9) 



For more details on the H-function as well as for a rich collection of its particular cases the reader is referred to [9] . 
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(n) 


X /d.o.t. 


10.7 


5.353 ± 1.431 


5.503 ± 0.171 


1.7/6 


13.8 


3.064 ± 0.298 


6.465 ± 0.151 


4.5/8 


13.9 


3.248 ± 0.250 


6.258 ± 0.112 


4.5/8 


16.6 


2.765 ± 0.087 


6.955 ± 0.071 


16.3/11 


18.2 


2.892 ± 0.551 


7 AO" ± 0.410 


2.1/8 


19.7 


2.942 ± 0.160 


7.603 ± 0.128 


7.8/11 


21.7 


2.57b ± 0.125 


7.838 ± 0.129 


15.4/12 


23.9 


o r"'70 1 ri oocT 


o.7o5 dz U.zz4 


13.5/11 


26.0 


2.725 ±0.113 


9.077 ± 0.114 


8.5/11 


27.6 


2.485 ± 0.101 


8.869 ± 0.136 


16.7/14 


30.4 


2.468 ± 0.150 


9.341 ± 0.261 


5.6/15 


38.8 


2.503 ± 0.076 


10.045 ±0.112 


11.9/14 


44. U 


z.o4y zt u. loo 


10.703 ± 0.223 


4.5/17 


52.6 


2.329 ±0.109 


11.457 ± 0.212 


10.4/19 


62.2 


2.392 ±0.105 


12.267 ±0.258 


17.5/18 


546.0 


1.840 ±0.034 


27.255 ± 0.292 


50.9/45 


TABLE I. 


Results of HNBD fits (WeibuU 


case, fc = 1) to the multiplicity distributions mea- 


sured in inelastic pp reactions [14]. The last 


row corresponds to the inelastic pp 


data of the UA5 


Collaboration. 






?7*-interval 


W (GeV) 




xVd-o.f. 


1< ??* < 2 


80 H- 115 


1.916 ±0.194 


1.6/13 




115 ± 150 


1.872 ± 0.085 


0.6/14 




150 ^ 185 


1.891 ±0.212 


1.5/14 




185 220 


1.841 ± 0.203 


0.8/14 


1< ?7* < 3 


80^ 115 


2.336 ± 0.166 


2.5/17 




115 -M 50 


2.116 ±0.143 


2.4/18 




150 ± 185 


2.041 ±0.143 


0.9/20 




185 ^ 220 


2.094 ±0.155 


0.6/21 


1< ??* < 4 


80 H- 115 


3.488 ± 0.385 


0.8/18 




115-^ 150 


3.100 ±0.180 


2.2/20 




150 ^ 185 


2.920 ± 0.226 


1.5/22 




185 ^ 220 


2.733 ±0.192 


1.6/22 


1< ?7* < 5 


80^ 115 


4.811 ± 0.392 


0.9/18 




115^ 150 


4.238 ± 0.277 


2.1/21 




150 ± 185 


4.011 ±0.288 


1.9/22 




185 ^ 220 


4.021 ± 0.334 


1.5/23 



TABLE II. Results of HNBD fits (Wcibull case, fc = 1) to the dccp-inclastic e+p multiplicity 
data of the HI Collaboration [15]. The average multiplicity (n) was fixed at its experimental value. 
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?7- interval 


a 


(n) 


xVd.o.f. 


Ar) = 3.0 


0.702 ± 0.020 


22.050 ± 0.546 


30.2/81 


Ar] = 5.0 


0.621 ± 0.018 


32.747 ± 0.780 


35.9/99 


full phase-space 


0.538 ± 0.014 


35.552 ± 0.834 


32.7/52 



TABLE III. Results of HNBD fits (log-normal case, p = 0) to the UA5 non-dilFractive pp mul- 
tiplicity data for |rj| > 2.5 at -^/s = 900 GeV [16]. The (p, ff) parametrization of the HNBD is 
discussed in the text. 



Experiment 

Aleph 

Delphi 

L3 

Opal 

Opal, 161 GeV 



0.202 ± 0.008 
0.201 ± 0.004 
0.210 ±0.009 
0.213 ±0.005 
0.235 ± 0.016 



21.094 ± 0.226 
21.353 ± 0.104 
20.742 ± 0.254 
21.349 ±0.103 
24.417 ± 0.423 



xVd.o.f. 
11.2/25 
32.4/24 
15.9/21 
17.2/25 
3.9/23 



TABLE IV. Results of HNBD fits (log-normal case, p = 0) to the e+e" 
plicity distributions at \/s = 91 GeV (top four rows) and 161 GeV [17]. 



full phase-space multi- 



J/- interval 


a 


(n> 


xVd-o.f. 


Ay = 0.5 


0.565 ± 0.046 


3.058 ±0.143 


3.5/18 


Ay = 1.0 


0.547 ±0.030 


6.426 ± 0.231 


9.3/32 


Ay = 1.5 


0.505 ± 0.018 


9.897 ±0.219 


26.3/38 


Ay = 2.0 


0.454 ± 0.021 


13.284 ± 0.293 


13.7/44 



TABLE V. Results of HNBD fits (log-normal case, p = 0) to the Aleph multiplicity data in 
central rapidity windows [18]. 




FIG. 1. HNBD fit (lognormal case, p = 0) to the non-difi^ractive pp multiplicity distribution in 
full phase-space at v'* = 900 GeV. 
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